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Abstract 


Objectives:  Develop  high  fidelity  acoustic  scattering  models  to  facilitate  the  detection,  local¬ 
ization,  and  characterization  of  military  munitions  found  in  ponds,  lakes,  rivers,  estuaries,  and 
coastal  ocean  areas.  Use  these  models  to  design  experiments,  interpret  collected  data,  identify 
features  that  can  be  used  in  training  classifiers. 

Technical  Approach:  Most  acoustic  scattering  models  are  well-equipped  to  compute  scat¬ 
tering  from  elastic  objects  in  free  space.  However,  unexploded  ordinance  (UXO)  reside  in 
complex  ocean  environments,  which  have  a  profound  effect  on  their  acoustic  response  to  sonar. 
Thus,  in  trying  to  model  the  response  of  UXO  to  sonar,  a  high-fidelity  model  must  account  for 
the  interaction  between  the  target  and  its  surrounding  environment.  While  the  standard  method 
of  solution  for  an  arbitrary  elastic  target  is  the  finite  element  method,  the  solution  of  the  scatter¬ 
ing  problem  in  the  surrounding  medium  is  best  handled  by  the  boundary  integral  equation,  as  it 
replaces  the  infinite  domain  problem  by  an  integral  over  the  surface  of  the  target.  Furthermore, 
the  boundary  integral  method  has  the  advantage  of  reducing  the  dimensionality  of  the  prob¬ 
lem  by  one.  In  contrast,  the  finite  element  method  is  not  well-suited  for  solving  the  scattering 
problem  in  the  surrounding  environment  due  to  difficulty  in  satisfying  the  radiation  condition. 
For  these  reasons,  we  solve  the  problem  of  scattering  from  an  elastic  target  in  a  complex  ocean 
environment  by  a  combination  of  finite  element  method  and  boundary  integral  method.  We 
use  the  finite  element  method  to  model  the  motion  of  the  target  by  essentially  computing  its 
impedance  matrix  in  vacuum,  and  the  boundary  integral  method  to  model  the  acoustic  field  in 
its  surrounding  medium.  The  two  solutions  are  coupled  by  satisfying  the  required  boundary 
conditions  on  the  surface  of  the  target.  This  results  in  a  model  that  treats  the  interaction  be¬ 
tween  the  target  and  its  surrounding  environment  exactly.  An  important  extension  of  this  work 
is  the  development  of  the  same  model  for  axially-symmetric  targets  with  substantial  speed  ad¬ 
vantages. 

Results:  During  the  past  three  years,  we  have  been  providing  benchmark-quality  solutions 
for  various  targets  to  researchers  within  SERDP  who  do  similar  type  of  modeling.  While  do¬ 
ing  this,  an  important  paid  of  our  work  has  been  to  validate  our  own  models  using  analytical 
solutions  when  available  and  other  well-tested  solutions.  We  validated  our  3D  and  axially- 
symmetric  models  in  free  space  using  the  analytic  solution  for  elastic  spheres  and  spherical 
shells.  We  also  validated  our  models  for  scattering  from  a  proud,  half-buried  and  a  fully -buried 
solid  sphere  by  comparing  our  results  with  those  of  the  T-matrix  method.  We  computed  the 
acoustic  color  (backscattered  target  strength  as  function  of  frequency  and  angle  of  incidence) 
for  the  aluminum  replica  of  a  UXO  (henceforth  aluminum  UXO),  the  Bullet- 105  and  the  How¬ 
itzer  shell  in  free  space  using  both  our  3D  and  axially-symmetric  versions  of  our  model.  We 
also  computed  the  acoustic  color  for  the  fully  proud  and  the  fully  buried  aluminum  UXO  and 
compared  our  results  with  measurements  and  those  produced  by  other  models.  Additionally, 
we  computed  the  acoustic  color  for  the  partially  buried  aluminum  UXO  and  again  compared 
our  results  with  other  finite  element  models  since  measured  results  for  this  case  was  not  yet 
available.  The  model  results  in  all  cases  agreed  with  each  other  and  and  with  the  measure¬ 
ments.  To  verify  that  our  model  indeed  accounts  for  multiple  scattering,  we  computed  the 
acoustic  color  for  a  tilted  cylinder  near  the  water-air  interface  and  showed  that  our  model  re¬ 
sults  were  in  total  agreement  with  measurements.  Finally,  to  show  that  our  3D  model  is  truly 
3D  and  can  handle  targets  of  arbitrary  shape,  we  used  it  to  compute  the  acoustic  color  for  the 
Bullet- 105  UXO  with  a  hole  drilled  on  its  side. 
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Benefits:  The  method  that  was  developed  has  several  advantages  over  currently-used  meth¬ 
ods.  The  most  important  ones  arc:  1)  The  method  is  inherently  broadband  since  the  stiff¬ 
ness  and  mass  matrices,  which  constitute  the  impedance  matrix,  are  independent  of  frequency. 
Therefore,  the  computation  of  these  matrices,  which  makes  up  the  most  numerically  intensive 
part  of  the  computation,  is  performed  once  for  all  frequencies.  2)  This  method  is  efficient 
because  it  requires  a  matrix  inversion  for  each  frequency,  but  not  each  angle  while  computing 
the  acoustic  color.  This  is  not  the  case  for  currently-used  methods,  which  must  solve  a  full 
finite  element  problem  for  each  frequency  and  each  angle  of  incidence.  3)  Since  this  method 
computes  the  target  impedance  matrix  in  vacuum,  the  same  impedance  matrix  can  be  used  in 
any  environment,  so  changing  the  environment  for  the  same  target  does  not  require  a  full  fi¬ 
nite  element  solution  of  the  problem.  4)  By  projecting  the  impedance  matrix  onto  the  surface 
nodes,  this  method  reduces  a  finite  element  problem  to  a  boundary  element  problem  with  far 
fewer  unknowns.  This  reduction  in  the  number  of  unknowns  enables  the  method  to  solve  a 
3D  problem  with  ease.  5)  It  provides  a  numerically  exact  solution  since  it  self-consistently 
couples  the  target  with  the  surrounding  environment.  6)  Due  its  modular  nature,  the  method 
easily  lends  itself  to  parallel  processing,  including  GPU  processing  and  the  application  of  the 
fast  multipole  technique. 

Objective 

The  objective  of  this  work  has  been  to  develop  high  fidelity  acoustic  scattering  models  that  can 
quickly  compute  the  acoustic  color  templates  (backscattered  target  strength  as  a  function  of  fre¬ 
quency  and  aspect  angle)  for  UXO  targets.  These  models  are  used  to  design  experiments,  interpret 
collected  data  and  particularly  identify  features  that  can  be  used  in  classification. 

Background 

Modeling  the  response  of  UXOs  to  a  sonar  signal  in  an  ocean  environment  belongs  to  a  large  class 
of  problems  referred  to  as  the  fluid-structure  interaction,  where  the  fluid  in  this  case  refers  to  the 
ocean  environment  and  the  structure  refers  to  the  UXO.  The  problem  of  determining  the  interac¬ 
tion  between  a  submerged  elastic  structure  and  its  surrounding  fluid  is  of  considerable  interest, 
particularly  in  underwater  acoustics  and  aeronautics  where  it  is  required  to  determine  the  acoustic 
field  about  an  arbitrary  three-dimensional  structure.  While  the  standard  method  of  solution  for  an 
arbitrary  elastic  structure  is  the  finite  element  method,  the  solution  of  the  reduced  wave  equation 
in  the  surrounding  medium  is  best  handled  by  the  boundary  integral  equation,  as  it  replaces  the 
infinite  domain  problem  by  an  integral  over  the  surface  of  the  submerged  structure.  Furthermore, 
the  boundary  integral  method  has  the  advantage  of  reducing  the  dimensionality  of  the  problem 
by  one.  In  contrast,  the  finite  element  method  is  not  well-suited  for  solving  the  wave  equation  in 
the  surrounding  fluid  environment  due  to  difficulty  in  satisfying  the  radiation  condition  as  well  as 
demands  on  the  mesh  size  and  the  difficulty  in  generating  the  fluid  mesh. 

For  these  reasons  the  problem  of  fluid- structure  interaction  is  perhaps  best  treated  by  a  combi¬ 
nation  of  finite  element  method  to  model  the  motion  of  the  structure  and  the  boundary  integral 
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method  to  model  the  acoustic  field  in  its  surrounding  medium,  where  the  coupling  between  the 
two  models  is  achieved  by  imposing  the  continuity  of  pressure  and  normal  particle  velocity  at  the 
surface  of  the  structure. 

The  coupled  finite  and  boundary  integral  method  has  been  used  by  several  authors  in  recent  years 
[l]-[  10]  and  (see  Amini  et  al.  [9]  and  the  references  therein).  The  main  differences  in  these  ap¬ 
proaches  are  the  particular  finite  element  package  and  boundary  integral  formulation  employed, 
the  numerical  approximation  used  and  the  details  of  the  method  of  coupling. 

In  this  work  [10],  we  use  the  coupled  finite  and  boundary  element  methods  to  compute  the  scattered 
acoustic  field  from  an  arbitrary  elastic  structure  in  an  arbitrary  medium,  characterized  by  a  Green’s 
function.  The  formulation  (also  known  as  the  impedance  formulation)  results  in  a  self-consistent, 
accurate,  and  numerically  efficient  solution  that  is  valid  in  any  environment.  In  modeling  the  struc¬ 
ture,  we  use  hexahedral  eight-node  elements  to  discretize  the  Galerkin  weighted  residual  equation 
and  in  modeling  the  acoustic  environment,  we  discretize  the  Helmholtz- Kirchhoff  equation  by 
expanding  pressure  and  normal  particle  velocity  in  piecewise  constant  basis  functions  over  quadri¬ 
lateral  elements  on  the  surface  of  the  structure.  The  finite  and  boundary  element  solutions  are 
coupled  by  imposing  the  continuity  of  pressure  and  normal  particle  velocity  on  the  surface  nodes 
of  the  structure.  An  important  part  of  this  work  is  the  development  of  the  coupled  finite  ele¬ 
ment/boundary  element  model  for  axially-symmetric  structures  with  non-axially-symmetric  load¬ 
ing,  which  reduces  the  3D  problem  to  several  uncoupled  2D  problems,  one  for  each  circumferential 
order  of  a  Fourier  expansion. 

We  refer  to  the  method  as  the  3D-IMP,  (IMP  for  impedance).  The  axially-symmetric  version 
of  the  model  is  referred  to  as  Axi-IMP. 

Methods 

In  this  section,  we  will  derive  the  3D  and  the  axially-symmetric  versions  of  the  IMP  method,  30- 
IMP  and  Axi-IMP.  We  will  validate  the  solution  for  a  sphere  in  free  space  and  two  half-space 
environments.  For  a  sphere  in  free  space  we  will  use  its  classical  partial  wave  solution  as  a  refer¬ 
ence  solution  and  for  a  sphere  in  two  half-space  environment,  we  will  use  the  T-Matrix  [11]  as  the 
reference  solution. 

The  fluid-structure  formulation 

We  use  the  fluid- structure  formulation,  based  on  the  work  of  Wilton,  Everstine,  Benthien  and 
Schenck  [3,  2,  1],  where  we  use  the  finite  element  technique  to  solve  the  Euler’s  equation  in  the 
elastic  target  and  couple  this  solution  with  the  solution  of  the  Helmholtz  equation  in  the  surround¬ 
ing  medium  by  using  pressure  loading  as  the  external  force  that  causes  displacements  in  the  elastic 
target.  A  more  complete  derivation  of  this  technique  is  provided  in  the  above  references,  but  to 
elucidate  the  discussions  that  follow,  we  summarize  the  derivation  below:  We  start  with  the  Euler 
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equation 


p^u  +  BTDBu  =  F,  (1) 

dtz 

Where  u  is  the  displacement  vector,  p  is  the  density,  F  is  the  external  force,  B  is  a  rectangular 
(6x3)  differential  operator  that  relates  the  strain  and  the  displacement  and  D  is  the  (6x6)  elasticity 
matrix.  The  displacement  in  the  elastic  target  can  be  expressed  in  the  finite  element  representation 
as 


N 


u(x)  =  YtuiVi(x), 

i=  1 


xeR, 


(2) 


where  R  is  the  region  the  object  occupies,  a,-  corresponds  to  the  components  of  the  displacement  at 
a  finite  number  of  points  in  the  object  called  nodes,  and  the  function  y/(-  forms  a  three-dimensional 
basis  for  piecewise  polynomial  interpolation  between  nodes.  The  finite  element  equations  can 
be  obtained  by  substituting  the  expansion  (2)  in  the  Euler  equation,  multiplying  it  with  t j/n  and 
integrating  over  the  volume  to  obtain 


d2U 

M~diT+KU=F' 


where 


Mmn  j  PWn  ^lr 


YlBTDBymdv, 


are  respectively  the  mass  and  stiffness  matrices.  If  the  time  dependence  is  harmonic  (e  1(01 )  with 
circular  frequency  ft),  the  above  equation  becomes 

(— ftrM  +  K)  U  =  F.  (3) 

On  the  ’wet’  surface  of  the  target,  the  Helmholtz-Kirchhoff  integral  can  be  written  as 

^p(x)  -  I p(xYJGj^ds'  =  -icop  I  v(x')G(x,x')ds'  +  Pinc(x),  x  e  5,  (4) 

where  p  is  the  pressure,  v  is  the  normal  particle  velocity  and  G(x,x')  is  the  environment  Green’s 
function.  Expanding  p  and  v  in  finite  element  basis  functions 

N  N 

p{x)  =  v(x)  = 

i=  1  ;=1 


and  substituting  them  in  (4)  gives 


where 


AP  =  BV  +  Pinc, 


8mn  f  d G (xn . xm )  ,  /  /  \  ,  / 

A mn  0  ds ,  Bmn  icop  /  G(xn,xm)as . 

2  JSn  dn  JSn 


(5) 

(6) 
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In  the  above  <j>,-  are  piecewise-constant  basis  functions  and  the  integration  is  performed  over  surface 
elements,  sn.  To  couple  the  elastic  and  the  fluid  equations,  we  note  that  at  the  surface  node  i  the 
force  is  due  to  pressure  loading 

Fi  =  -  J  p(q) nq  ■  \ifj(q)dSq. 

On  the  other  hand 


p(q)  =  X>A(<?)- 

7=1 


From  these  two  equations,  we  obtain 

F  =  -I 7p, 


where 

Ljj  —  J (j)i(q)  yfj(q)  ■  nqdSq. 

The  continuity  of  normal  displacement  requires  that  the  normal  velocity  in  the  fluid  should  be 
related  to  the  normal  displacement  in  the  structure  as 


V  =  —icoLu. 


Substituting  this  expression  in  (3)  and  (5),  we  find 

A  P  =  -i(0BLu  +  Pinc, 

(— «2M  +  K)u  =  — I/P. 

Eliminating  u  in  the  above  two  equations  gives  the  pressure 

P=  (A-zc»BL(-®2M  +  K)_1i/)  '  Pinc,  (7) 

and  normal  particle  velocity  on  the  surface 

V  =  /®L(— ©2M  +  K)_1l/P  (8) 

Equation  (7)  relates  the  total  pressure  field  on  the  surface  of  the  object  to  the  incident  field  and 
a  matrix  containing  information  about  the  surrounding  environment  in  the  matrices  A  and  B  and 
properties  of  the  elastic  object  in  matrices  L,M,  and  K.  The  latter  three  matrices  are  computed 
using  a  finite  element  model  and  the  former  two  and  pmc  are  computed  using  a  propagation  model. 
Once  P  is  computed  using  (7),  (8)  is  used  to  compute  the  normal  velocity,  V.  Using  the  surface 
pressure  and  normal  velocity,  the  scattered  field,  P\  is  computed  at  any  arbitrary  point  using  the 
Helmholtz-Kirchhoff  integral 

Ps  =  aTP-\-bTV,  where  an  =  f  ,x Ids',  bn  —  icop  f  G(xnx')dsr.  (9) 
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The  most  numerically  intensive  part  of  computing  (7)  is  computing  the  inverse  of  the  impedance 
matrix,  (— ctrM  +  K),  since  the  dimensions  of  M  and  K  are  equal  to  the  number  of  volume  ele¬ 
ments,  which  includes  interior  as  well  as  surface  elements.  What  proved  to  be  a  crucial  step  in  this 
method  is  the  use  of  modal  decomposition  to  invert  the  impedance  matrix,  based  on  the  work  of 
Benthien  and  Schenck  [3].  Introducing  the  matrix  E  whose  columns  are  the  normal  modes  of  the 
structure, 


KE  -  MEQ, 

where  Q.  is  the  diagonal  matrix  of  in  vacuo  eigenfrequencies  of  the  structure.  Since  the  modes  can 
be  normalized  so  that  it  is  M-orthogonal  and  K-normalized,  we  can  write 

ErME  =  I. 


Then  it  can  be  easily  verified  that 

(-gtM  +  K)”1  =E(-arI  +  n)-1Er.  (10) 

Thus,  the  impedance  matrix  can  be  inverted  for  all  frequencies,  except  those  that  are  the  in  vacuo 
eigenfrequencies  of  the  structure.  If  a  desired  frequency  happens  to  coincide  with  one  of  the 
eigenfrequencies,  linear  interpolation  between  neighboring  frequencies  can  be  used.  Through  nu¬ 
merical  experimentation  we  have  found  that  the  above  result  is  accurate  if  the  maximum  operating 
frequency  is  lower  than  the  eigenfrequency  of  the  highest  mode.  This  can  simply  be  accomplished 
by  specifying  the  number  of  modes  to  be  computed  by  the  finite  element  model  and  making  sure 
that  the  eigenfrequency  of  the  highest  mode  is  larger  than  the  highest  frequency  of  interest.  The  use 
of  (7),  (8),  (10)  in  (5)  completely  solves  the  scattering  problem  in  an  environment  whose  Green’s 
function  is  specified  in  (6). 

Modification  for  axially-symmetric  objects  in  free  space 

By  a  process  that  will  be  described  below,  it  can  be  shown  that  the  solution  of  scattering  from  an 
axially-symmetric  object  with  non  axially-symmetric  loading  is  equivalent  to  solving  a  3D  problem 
by  a  series  of  2D  problems,  which  results  in  significant  reduction  in  numerical  cost.  Referring  to 
Fig.  (1)  we  note  that  if  an  axially-symmetric  object  is  insonified  by  an  incident  field  that  is  along 
its  axis  of  symmetry,  the  field  on  a  line  such  as  A,  which  is  formed  by  the  surface  of  the  object  and 
a  plane  perpendicular  to  the  axis  of  symmetry,  is  constant. 
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Figure  1:  The  geometry  of  an  axially-symmetric  object. 

However,  if  the  incident  field  is  not  along  the  axis  of  symmetry  of  the  object,  the  field  along  A  is 
periodic  in  azimuthal  angle  (j).  Since  this  field  is  only  a  function  of  (j).  which  can  only  happen  for 
an  axially-symmetric  object  in  which  case  A  is  a  circle,  it  can  be  expanded  in  an  azimuthal  Fourier 
series 


P  (x)  =  £p”(r,z)cos(n0)+p"(r,z)sin(n0),  (11) 

n= 0 

Since  the  Fourier  coefficients,  pn,  are  not  a  function  of  (j) ,  they  represent  an  axially-symmetric  load 
(incident  field).  Therefore,  the  above  expansion  is  a  sum  of  axially-symmetric  solutions  when  the 
loading  is  not  axially-symmetric.  If  the  loading  is  axially-symmetric,  then  only  the  term  n  =  0 
contributes.  Substituting  the  above  Fourier  expansion  for  the  pressure  and  normal  particle  velocity 
into  (4)  and  switching  the  order  of  sums  and  integrals  gives 

-t  oo 

2  L  [p"(r,z)cos(n0)+jp"(r,z)sin(n0)] 

^  n=0 


~Y  [  [Pne  (rV)  cos(»0')  +Po  sin(n0')l  dG[X'X  ^  ds' 

n=0J 


dn' 


(12) 


=  -icop  Y,  /  [ve  ( r/,z/)c°s(77^ >')  +v"  (r',z')  sin G(x,x')ds'  +  pinc(x). 

n=—oo  J 

Applying  the  integral  operator 


c2k 


cos (£(j))  d(j), 
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to  the  above  equation  results  in  the  boundary  element  integral  equation  for  the  even  azimuthal 
coefficients  given  by 


\pl(r,z)  -  I ple(r',z  )hle(r:z-rl:z  )r'dr' 

—  —i(op  f  vee  (/,  z')gl(r,  z;  r  ,z)r  dF'  +  Pinc e}  ’  ^  . 


(13) 


where 


PincM )  =  J  Pinc(X)  COS(£0)  d(j>, 


and 

1  r2  71  p2k 

g[  (r,z;r',z)  =  —  j ^  G(x,x')  (cosf(^/-<^)+cosf(^/  +  0)) 

(14) 

1  rlTZ  rlK  Xf) 

hi(r,z-,r,z)  =  —  ^  ^  (cosf(0/-<^)+cosf(^/  +  <^))  ##', 

where  e  =  2  for  t  =  0  and  1  otherwise.  The  boundary  element  equation  for  the  odd  azimuthal 
coefficients  is  obtained  by  applying  the  operator 

/•2k 

/  sin (£<j>)  d<\ >, 

Jo 


to  (12).  It  is  given  by 


where 


^Po(pz)  —  j p£0(r  ,z)he0(r,z-,r  ,z!)r'dT' 

=  -i(Dp  j vC0(r,z)gi{r,z;r',z)r'dr'  +  I  inc°^ 


1  r  £71  r£7 1 

g£0(r,z;r,z)  =  —  G(x,x')  (cos£(</  -  <p)  -  cos£(<j>'  +  0))  d<f>  dty' } 

[cosf(^/  —  (j))  —  cos £($'  +  $))  d(j)  d(j)f , 


1  P  f2’ <>G(x,x') 


ho(^y,z)  =  -h  h  M 


and 


(15) 


(16) 


/■2tt 

Pinc0(x)=  Pmc(x)sin(^)d0. 
./o 
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If  G  is  a  function  of  the  difference  of  the  angles  ((/)'  —  (j)),  as  is  the  case  for  the  free  space  Green’s 
function,  (14)  and  (16)  reduce  to 


[•2n 

ie  (r,z-J,z!)  =ge0  (r,r,r',z)  =  J  G(x,x/)cos(^/)#', 


hi  (r,z;r',z)  =  he0  (r,z;r,z)  =  J  ~  cos  (ttf)dtf , 


(17) 


However,  for  a  general  Green’s  function  these  equations  will  be  coupled  in  azimuthal  order  and 
require  n2  evaluations  of  the  double  angular  integrals  in  (14)  and  (16).  In  this  case,  the  axially 
symmetric  boundary  element  formulation  has  no  computational  advantage  over  its  3D  counterpart. 
Discretizing  the  above  equations  as  was  done  before  and  assuming  piecewise  constant  fields,  we 
get  for  each  azimuthal  order  i 


A  p  -K  v  +pinci 


where  the  matrix  elements  are  given  by 


A1  )  —  — 

nmn  /  —  o  umn 

e.o  Z 


j  h^o(rm,zm-yn,z!nyndT' , 


(18) 


(19) 


=  -icop 


(20) 


where  all  variables  (r,z;  r' .z')  in  (19)  and  (20)  are  evaluated  on  the  surface.  Note  that  (18),  (19)  and 
(20)  are  the  axially  symmetric  counterparts  of  (5)  and  (6).  These  equations  are  valid  in  the  fluid 
and  on  the  wet  surface  of  the  structure.  The  axially  symmetric  counterpart  of  (3)  for  the  structure 
is  given  by 

(-w2Mf  +  Kfju'-Ff,  (21) 

where 

(ML)eo  =  fp  vl  K1K1  ¥„  rVrVf, 

and 


Win 


\Cf°]BT DB\Cp°]  \ffn  r'dT'dtf. 


The  displacement  basis  functions  are  given  by 


’  Wr  ’ 

oo 

’  Yr  ’ 

e 

’  Yr  ’ 

Vm  = 

Yz 

=  £  [Cf] 

Yz 

+  [cfi 

Yz 

.  Ye  _ 

£=0 

m 

.  Ye  _ 

m ,1 

.  Ye  _ 
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and 


cos  (£<j>)  0  0  sin  (£(j>)  0  0 

[Cf]  =  0  cos  (£(/))  0  ,  \C°A  =  0  sin  (£<j>)  0 

0  0  —  sin(70)  0  0  cos(£(/)) 

The  combination  of  (18)  and  (21)  results  in  similar  equations  as  (7)  and  (8)  for  the  axially  sym¬ 
metric  case 


Pe,o  =  Ko  ~  ^BloL  [-CO2  M{0  +  K{0 


Lr  Pfnce^ 


Ve,0  =  ia>  L  (-co2  Mji0  +  <0  )  L%. 


Note  that  for  an  axially- symmetric  problem  the  sizes  of  the  matrices  Me  o  and  Kre  o  are  small 
enough  that  the  matrix  inversion 


-(02Mio  +  Kio)  , 


can  be  performed  directly  without  the  use  of  the  modal  expansion  described  by  (10). 


The  integrals  in  (19)  and  (20)  are  performed  over  the  generating  curve  F.  If  F  is  divided  into 
N  linear  elements,  an  integral  of  the  form  given  in  (19)  and  (20)  over  a  single  element  n  can  be 
written  as 


7/77/7  - 


+  fe(rm,Zm',r,z)rJ  1  +  (^j  dz, 


where  the  arc  length  dT  —  \J  1  +  ( dr/dz )2.  Transforming  the  above  integral  to  a  local  coordinate 
system  described  by 

-  fJn+\  +r'n  ,  <+l 

-  ^n+ 1  <-n  ^/z+l  Zn  £ 

results  in 

1  1 

7/77/7  =  ^  J  Qmn  ( rm 7 Zm !  )d^n,  (24) 

where 

Qe,nn  (rm,Zm\Zn)  -  fe(rm,Zm\ £n)  {rn+1  +  r'n  +  (rn+1  -  r'n)  Q  \J (rJn+l-r,n)2  +  (z,n+l-Z,n)2. 
Using  two-point  Gauss  quadrature,  the  integral  in  (24)  can  be  evaluated 

lynn  =  Qmn  (j miZm'i  ~~  1/V^3^  +  Qmn  Z.m\  1/V^3^  • 

The  angular  integrals  in  (17)  are  evaluated  using  multi-point  Gaussian  quadrature. 
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Validation 


Due  to  its  well-known  analytic  solution,  the  sphere  provides  the  best  method  for  validating  a  nu¬ 
merical  solution.  Here  we  use  a  solid  steel  sphere  and  an  aluminum  spherical  shell  to  benchmark 
Axi-IMP.  Figure  (2)  shows  a  comparison  between  the  backscattered  pressure  amplitude  from  a  0.5 
m  solid  steel  sphere  computed  using  Axi-IMP  and  the  classical  partial  wave  solution  as  a  function 
of  frequency. 


Figure  2:  The  scattering  amplitude  as  a  function  of  frequency  for  a  0.5  m  solid  steel  sphere.  The 
solid  line  represents  the  partial  wave  solution  and  the  red  dash-dot  line  represents  the  coupled 
finite-boundary  element  solution. 

A  comparison  of  the  backscattered  scattering  amplitude  computed  using  Axi-IMP  and  the  partial 
wave  solution  for  an  air-filled  aluminum  spherical  shell  is  shown  in  Fig.  (3).  The  radius  of  the 
sphere  in  this  case  is  0.5  m  and  the  shell  thickness  is  5  cm.  In  both  cases,  there  is  excellent 
agreement  between  our  solutions  and  the  benchmark  solutions. 


01  23456789  10 

Frequency  (kHz) 


Figure  3:  The  scattering  amplitude  as  a  function  of  frequency  for  a  0.5  m  air-filled,  aluminum 
spherical  shell.  The  solid  line  represents  the  partial  wave  solution  and  the  red  dash-dot  line  rep¬ 
resents  the  coupled  finite-boundary  element  solution. 
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We  validated  the  3D  version  of  our  model,  3D-IMP,  by  using  spheres  in  free  space  and  were  able 
to  obtain  the  same  excellent  agreement  with  the  partial  wave  solutions  as  are  shown  in  Figs.  (2) 
and  (3).  The  Green’s  function  used  in  the  above  two  cases  is  the  free  space  Green’s  function  given 
by 


Jky/  ( x-x')2+{y-y,)2+{z-z ')2 

G  (x; x')  = - .  = .  (25) 

47 Vy/  (x-x')2  +  (y  —  y')2  +  (z  —  z')2 

Another  class  of  problems  with  a  slightly  more  complicated  Green’s  function,  which  also  lends 
itself  to  an  analytical  solution  is  the  case  of  an  acoustically  soft  or  rigid  hemisphere  on  an  infinite 
plane  satisfying  the  same  boundary  conditions.  This  is  also  referred  to  as  an  infinite  plane  with  a 
hemispherical  boss.  The  analytical  solution  is  obtained  from  the  partial  wave  solution  using  the 
method  of  images.  The  Green’s  function  used  in  the  3D-IMP  formulation  for  the  acoustically  soft 
case  with  the  infinite  plane  at  z  —  0  is 

GD  (x;x')  =  G\  (x;x')  -  G2  (x;x') ,  (26) 


where 


pik\J  (x-x')2+(y-y')2+(z-z')2 

G\  (x;x')  = - e 

47T  \J  (x  —  x')2  +  {y  —  y')2  +  (z  —  z')2 

Jk\j  {x—x')2+{y—y')2jr(z+z')1 

G2  (x;x')  = - * 

4k\J{x-x')2  +  {y-y')2  +  (z  +  z')2 

For  the  rigid  case  the  Green’s  function  is  given  by 

Gn  (x;  x')  =  G\  (x;  x')  +  G2  (x;  x') .  (27) 


These  Green’s  functions  guarantee  that  the  total  field  vanishes  on  the  infinite  plane  for  the  acous¬ 
tically  soft  case  and  its  normal  derivative  vanishes  on  the  same  surface  for  the  acoustically  hard 
case.  The  boundary  element  method  then  enforces  the  same  boundary  conditions  on  the  surface  of 
the  hemispherical  boss,  resulting  in  a  self-consistent  solution.  In  this  example,  the  hemispherical 
boss  has  a  0.5  m  radius  and  the  problem  is  solved  for  a  plane  wave  making  a  30°  angle  with  the 
positive  jc-axis  for  two  frequencies  of  5  kHz  and  10  kHz.  The  scattered  field  is  computed  as  a  func¬ 
tion  of  the  receiver  angle.  The  surface  of  the  hemispherical  boss  was  meshed  by  approximately 
8000  quadrilateral  elements.  The  solutions  are  shown  in  Fig.  (4),  where  they  are  compared  with 
the  corresponding  exact  partial  wave  solutions. 
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Figure  4:  The  scattering  amplitude  as  a  function  of  the  receiver  angle  for  a  hemispherical  boss. 
The  scattering  geometry  is  shown  in  the  inset  in  the  top  left  panel.  The  incident  plane  wave  makes 
an  angle  y  =  30°  with  the  x-axis  in  all  cases.  The  receiver  angle,  d  varies  from  zero  to  180°.  The 
top  and  bottom  panels  show  the  results  for  the  acoustically  soft  and  rigid  cases,  respectively.  The 
left  and  right  panels  are  for  the  incident  frequency  of  5  kHz  and  10  kHz.  The  solid  lines  are  for  the 
exact  partial  wave  solutions  and  the  broken  red  lines  are  for  coupled  3D-1MP  solution. 


The  next  order  of  complexity  arises  when  the  ocean  environment  is  composed  of  two  fluid  half 
spaces:  a  water  half  space  over  a  fluid  bottom  half  space.  In  most  practical  cases  the  ocean  environ¬ 
ment  is  modeled  this  way,  so  it  is  appropriate  to  validate  this  solution  for  this  type  of  environment. 
Since  the  T-Matrix  solution  is  designed  for  spherical  targets,  in  this  validation  process  we  will  use 
a  solid  sphere.  We  will  consider  three  cases,  fully  proud,  fully  buried  and  partially  buried.  As  will 
be  seen,  the  latter  case  proves  to  be  quite  challenging. 

If  the  waveguide  is  composed  of  a  water  half-space  with  sound  speed  c \  and  density  p i  over  a 
bottom  half-space  with  sound  speed  C2,  density  P2,  and  interface  at  z  —  0,  the  Green’s  function  is 
given  by 


G(r,z;rs,zs) 


G{knz;zs)Jo(kr\r 


-rs\)krdkr, 


(28) 
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where  rs  and  zs  are  the  source  range  and  depth  and  r  and  z  are  the  receiver  range  and  depth,  Jo  is 
the  Bessel  function  of  order  zero  and  the  spectral  Green’s  function  G(kr,z',zs )  for  a  source  located 
in  the  water  is  given  by 


' 

G{kr,z;zs )  =  < 


iReikz  i  (z+z0  ieikz  i  I 

4  nkZl  +  47 vkZl 

iTeikziZse~ikz2z 


AnkZl 


for  z  >  0, 

for  z  <  0, 


(29) 


where  fcZl  —  \Jk\  —  k},  kZ2  =  k\  —  kj,  ft)  is  the  circular  frequency  and  the  reflection  and  transmis¬ 
sion  coefficients  are  given  by 

kZ\  Pi  ~  kZ2p\  2kz !  P2 


/?  = 


kZAp2  +  kZ2p\ ' 


T  = 


ft) 


7.7)  )  i  1)2- 

%P2+%Pl  ci 


The  above  Green’s  functions  are  used  to  compute  the  incident  field,  the  matrices  A  and  B,  given 
in  (6)  and  the  scattered  field  given  in  (9).  The  matrices  A  and  B  represent  the  interaction  between 
surface  elements,  where  one  element  plays  the  role  of  the  source  and  the  other  the  role  of  the 
receiver  in  (29)  and  (30).  Similarly,  in  computing  the  scattered  field,  the  surface  elements  play  the 
role  of  the  sources.  For  this  reason,  if  the  entire  target  or  any  part  of  it  is  in  the  bottom  half  space, 
we  need  a  spectral  Green’s  function  for  a  source  in  the  bottom,  which  is  given  by 

iTeikzize-ikz2z° 

—  for  z>  0, 


Ajik 


G(kr,z;zs)  =  { 


Z2 


—iRe~ikz  2(Z+Zj)  ieikz2\z-zs\ 

-t - — -  for  z<0, 


4  nk7 


4  nkv. 


vZ2  'tA'Z2 

where  now  the  reflection  and  transmission  coefficients  are  given  by 

^ziP2~kZ2pi  _  2  kZ2p\ 


R  = 


kZlp2  +  kZ2pi  kZlp2  +  kZlpi 

By  substituting  (29)  and  (30)  in  (28),  we  obtain  for  the  source  in  the  water 

iReik^z+zP 


Ank7 


-J0(kr\r-rs\)krdkr  +  gw(ki,r,z;rs,Zs)  for  z  >  0, 


G{r,z;rs,zs)  =  { 


iTeikz  iz°e~ikz2z 


/  o 


47 ik7 


-J0(kr\r -  rs\)k,-dkr 


and  for  the  source  in  the  bottom 


iTeikzPe~ikz2Zs 


-Jo(kr\r  —  rs\)krdkr 


for  z  <  0, 


for  z  >  0, 


G(r,z;rs,zs)  =  { 


Z2 


K  JO 


-iRe~ikz  2(z+Zi) 
Anky^ 


Jo (kr\r  —  rs | )krdkr  +  gw (k2.J  r, z\  rs . zs)  for  z  <  0, 


(30) 


(31) 


(32) 
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where 


gw(k,r,z;rs,zs)  = 


eik\/  ( r-rs)2+(z-zs )2 

An  \/(r  -  rs)2  +  (z  -  z.*)2  ’ 


is  the  free  space  Green’s  function,  which  represents  the  direct  wave1  and  can  also  be  represented 
by  the  Sommerfeld-Wyle  integral  in  the  spectral  domain 


gw(k,r,z;rs,zs )  =  ^  J 


oo  Jkz\z-Zs\ 


-JQ{kr\r  -  rs)krdkr. 


The  computation  of  matrices  A  and  B  requires  n2  evaluations  of  the  Green’s  functions,  where 
n  is  the  number  of  surface  elements.  For  a  typical  target,  we  have  anywhere  between  6  to  20 
thousand  surface  elements,  which  means  that  the  computation  of  these  matrices  requires  on  the 
order  of  one  hundred  million  function  evaluations.  This  becomes  a  daunting  numerical  task  when 
the  Green’s  functions  are  represented  in  terms  of  spectral  integrals  as  in  (31)  and  (32).  For  a 
target  located  entirely  in  the  water,  the  Green’s  function  is  given  by  the  top  equation  in  (31)  and 
for  one  located  entirely  in  the  bottom  it  is  given  by  the  bottom  equation  in  (32).  We  note  that 
in  these  two  cases  the  Green’s  functions  are  a  function  of  the  sum  of  source  and  receiver  depths, 
(z  +  Zj),  and  the  difference  of  their  ranges,  \r  —  rs\.  So  in  these  two  cases,  the  Green’s  functions 
are  a  function  of  two  parameters  as  opposed  to  four  (r,z',rs,zs),  which  allows  for  the  use  of  2D 
interpolation.  For  the  cases  of  fully  proud  and  fully  buried  targets,  we  computed  the  Green’s 
functions  on  a  200  x  200  grid  as  a  function  of  the  sum  of  source  and  receiver  depths  and  the 
difference  of  their  ranges  around  the  target,  and  computed  their  values  at  the  surface  elements  using 
two-parameter  interpolation.  This  essentially  reduced  the  number  of  function  evaluations  from 
108  to  104  and  proved  to  be  a  crucial  step  in  our  ability  to  compute  the  acoustic  color  for  targets 
of  interest.  Figure  (5)  shows  the  interaction  between  surface  elements  and  their  corresponding 
Green’s  functions.  For  the  case  of  a  partially  buried  target,  the  situation  is  more  complicated. 
Referring  to  Fig.  (5),  in  this  case  the  Green’s  functions  are  a  function  of  source  depth,  zs,  receiver 
depth,  z,  and  the  difference  between  their  ranges,  | r  —  rv | ,  which  would  require  a  3-parameter 
interpolation.  Three-parameter  interpolation  is  also  numerically  expensive  and  does  not  provide 
a  computational  advantage.  To  deal  with  this  situation,  we  divide  the  interaction  matrix  (A  or  B) 
into  four  sub-matrices,  depending  on  the  locations  of  sources  and  receivers.  This  is  shown  by  the 
cartoon  in  Fig.  (6).  We  compute  sub-matrices  1  and  4  by  using  interpolation  as  we  did  for  the 
cases  of  fully  proud  and  fully  buried  targets,  respectively.  This  is  also  shown  in  the  top  two  panels 
of  Fig.  (5).  We  compute  sub-matrix  2  directly  and  sub-matrix  3  from  sub-matrix  2  by  applying  a 
special  form  of  reciprocity  according  to  the  equations 


83(kr,z,zs )  =  —g2{knz:zs)T: 

Pi  T 

dg3{kr,z,zs)  _  _pi  / dgi{knz,zs)  \ 

dx  P2  V  dx  J  ’ 

dg3(kr,Z.,Zs)  __plf dgl{kr,Z:Zs)\T 

dy  Pi  V  dy  ) 

dg3{kr,z,zs)  =  _.k  pi  ( dg2{kr,z,zs)\ 
dz  72  Pi  V  dz  ) 


T 


'The  direct  wave  is  a  wave  that  does  not  interact  with  the  boundaries  of  the  waveguide. 
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Figure  5:  The  interaction  between  surface  elements  and  their  corresponding  Green ’s  functions  for 
a  fully  proud,  fully  buried  and  partially  buried  target. 


where  g2  is  the  spectral  Green’s  function  for  sub-matrix  2,  g^  is  the  same  for  sub-matrix  3  and  the 
superscript  T  indicates  matrix  transpose.  So  in  doing  this,  we  have  reduced  the  computation  from 
()((£  +  m)2)  to  approximately  O(mT),  where  i  x  m  is  the  size  of  sub-matrix  2.  This  gave  us  the 
ability  to  compute  scattering  from  a  partially  buried  target  in  a  reasonable  time. 
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Figure  6:  The  interaction  matrix  is  divided  into  four  sub-matrices  depending  on  the  locations  of 
sources  and  receivers. 


To  validate  our  solution,  we  computed  scattering  from  a  fully  proud,  fully  buried  and  partially 
buried  solid  sphere  and  compared  our  results  to  those  obtained  by  the  T-Matrix  method.  The  T- 
Matrix  method  is  a  semi-analytic  formulation  [12],  based  on  the  analytic  solution  of  the  wave 
equation  for  a  sphere.  For  this  reason,  it  produces  almost  exact  solutions  for  spherical  objects.  We 
used  a  version  of  the  T-Matrix  method  that  is  designed  to  compute  scattering  from  a  sphere  in  a 
multi-layer  environment  [11].  In  this  simulation,  the  sphere  is  made  of  tungsten  carbide  with  a 
radius  of  6.35  mm.  The  environment  consisted  of  a  water  half  space  on  top  of  a  fluid  bottom  half 
space.  The  water  sound  speed  is  1485.4  m/s  with  a  density  of  1000  kg/m3,  and  the  bottom  sound 
speed  and  density  are  1655  m/s  and  1890  kg/m3,  respectively.  We  computed  the  backscattered 
scattering  amplitude  as  a  function  of  ka  for  a  proud,  buried  and  half  buried  sphere  at  a  backscattered 
angle  of  50°,  measured  in  the  counterclockwise  direction  from  the  interface,  where  k  =  (o/c\  and 
a  is  the  radius  of  the  sphere.  The  scattering  geometry  for  each  case  is  shown  in  the  cartoon  in  the 
left  panel  of  each  figure,  where  the  arrows  only  indicate  the  angle  of  incidence  and  scattering,  but 
do  not  represent  the  ray  paths.  This  is  an  exact  model  and  includes  all  multi-paths.  In  Fig.  (7)  we 
show  a  comparison  of  our  result  with  that  of  the  T-Matrix  solution  for  a  fully  proud  sphere. 
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Figure  7:  Comparison  of  backscattered  scattering  amplitude  as  a  function  ofkafor  a  solid,  proud 
tungsten  carbide  sphere  computed  using  the  T-Mcitrix  method  (black  line)  and  the  3D-IMP  method. 
The  scattering  geometry  is  shown  in  the  figure  on  the  left.  The  incident  field  is  a  plane  wave  and 
0  =  50°.  The  receiver  is  at  one  meter  from  the  center  of  the  sphere  in  the  direction  of  the  incident 
field. 


The  scattering  amplitude  as  a  function  of  ka  for  the  same  sphere  buried  in  the  bottom  such  that  its 
center  is  8  mm  below  the  interface  is  compared  with  the  T-Matrix  solution  in  Fig.  (8). 
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Figure  8:  Comparison  of  backscattered  scattering  amplitude  as  a  function  ofkafor  a  buried  tung¬ 
sten  carbide  sphere  computed  using  the  T-Matrix  method  (black  line)  and  the  3D-IMP  method. 
The  scattering  geometry  is  shown  in  the  figure  on  the  left,  where  0  =  50°  and  the  receiver  is  at  one 
meter  from  the  center  of  the  sphere  in  the  direction  of  the  incident  field. 


The  scattering  amplitude  as  a  function  of  ka  for  the  same  sphere,  half  buried  is  compared  with  the 
T-Matrix  solution  in  Fig.  (9). 
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Figure  9:  Comparison  of  backscattered  scattering  amplitude  as  a  function  ofkafor  a  half  buried 
tungsten  carbide  sphere  computed  using  the  T-Matrix  method  ( black  line)  and  the  3D-IMP  method. 
The  scattering  geometry  is  shown  in  the  figure  on  the  left,  where  0  =  50°  and  the  receiver  is  at  one 
meter  from  the  center  of  the  sphere  in  the  direction  of  the  incident  field. 

The  agreement  between  the  two  models  in  all  the  above  cases  is  excellent.  The  minor  differences 
between  the  two  models,  particularly  in  the  case  of  the  half-buried  sphere,  are  hard  to  explain 
since  both  models  are  susceptible  to  numerical  error  and  the  partially  buried  case  is  also  more 
challenging  to  the  T-Matrix  method.  Errors  at  higher  frequencies  could  also  be  because  we  use 
fewer  elements  than  necessary.  Increasing  the  number  elements  increases  computational  time  and 
dealing  with  it  is  an  issue  that  we  plan  to  address  in  the  future. 

The  advantages  of  the  IMP  method 

Most  problems  of  interest  to  SERP  are  numerically  too  demanding  to  be  solved  on  a  regular  desk¬ 
top  computer  by  conventional  finite  element  model  that  treats  the  target  as  3D  and  includes  the 
environment  as  part  of  the  computational  domain.  The  development  of  the  model  AxiScat  [13] 
was  an  effort  to  reduce  computational  complexity.  This  model,  which  is  in  use  by  some  inves¬ 
tigators  working  for  SERDP  computes  scattering  from  an  axially- symmetric  target  in  free  space 
for  a  given  incident  plane  wave.  Its  limitations  are  that  the  target  must  be  axially-symmetric  and 
the  environment  that  embeds  it  must  be  free  space.  But  in  most  simulations  involving  UXOs, 
the  environment  is  modeled  by  two  half-spaces.  The  presence  of  the  interface  causes  singly  and 
multiply-reflected  rays  to  also  ensonify  the  target.  To  approximate  the  solution  in  this  environment 
using  AxiScat,  scattering  due  to  these  rays  is  computed  individually  and  added  to  the  scattering  that 
is  caused  by  the  direct  ray.  This  is  an  approximate  solution  since  in  solving  the  scattering  problem, 
it  is  always  assumed  that  the  target  is  in  free  space  and  the  effects  of  the  interface  is  added  later, 
which  ignores  coupling  effects.  This  solution  is  referred  to  as  the  Hybrid  solution  since  it  uses  a 
combination  of  finite  elements  and  ray  theory  to  solve  the  problem.  In  contrast,  our  solution  is 
numerically  exact  since  the  environment  is  integrated  into  the  solution  through  the  environment 
Green’s  function  and  it  can  be  any  environment.  For  this  reason,  part  of  our  collaborative  work 
with  other  SERDP  researches  has  been  to  provide  benchmark  solutions.  Other  advantages  of  the 
method  are  listed  below: 
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1 .  The  most  important  advantage  of  the  IMP  technique  is  its  speed  and  efficiency  in  computing 
the  scattered  field  for  a  range  of  frequencies.  While  the  Hybrid  method  solves  the  scattering 
problem  for  every  frequency,  the  IMP  method  uses  the  finite  element  technique  to  compute 
the  matrices  C,  M,  K,  D  and  E  in  vacuum  independent  of  frequency  and  uses  (7)  to  com¬ 
pute  the  scattered  field  for  any  desired  frequency.  This  technique  is  much  more  efficient  in 
computing  the  acoustic  color. 

2.  The  other  advantage  of  the  IMP  technique,  which  is  enormously  useful  in  computing  acous¬ 
tic  color  plots,  is  that  it  only  requires  matrix  inversion  for  every  frequency,  but  not  for  every 
angle.  This  is  because  in  (7)  the  incident  field,  which  contains  the  angle  information  appears 
on  the  right-hand  side.  In  contrast,  since  the  Hybrid  method  solves  a  finite  element  problem 
for  every  frequency  and  every  incident  angle,  it  requires  matrix  inversion  at  every  frequency 
and  angle. 

3.  The  IMP  technique  completely  decouples  the  surrounding  environment,  represented  by  ma¬ 
trices  A  and  B  from  the  elastic  target  represented  by  matrices,  C,  M,  K,  D  and  E.  Once  these 
matrices  are  computed  in  vacuum,  they  can  be  used  in  any  environment.  One  can  therefore 
construct  a  library  of  impedance  matrices  for  all  targets  of  interest  and  compute  scattering 
from  them  in  any  desired  environment.  Furthermore,  by  solving  the  finite  element  prob¬ 
lem  in  vacuum,  there  is  no  need  to  use  PMLs  or  similar  techniques  to  impose  the  radiation 
condition.  These  techniques  are  often  tricky  and  seldom  unreliable. 

4.  The  IMP  technique  self-consistently  couples  the  target  with  the  surrounding  environment 
and  provides  a  numerically  exact  solution  of  the  problem  regardless  of  whether  the  target 
is  in  free  space,  a  multi-layered  medium  or  a  waveguide.  Most  conventional  finite  element 
techniques,  on  the  other  hand,  compute  scattering  in  a  small  computational  domain  around 
the  target  and  can  only  treat  the  presence  of  boundaries  outside  of  this  domain  using  a  single¬ 
scatter  approximation,  which  ignores  multiple  scattering.  As  is  shown  in  the  next  section, 
this  solution  becomes  highly  inaccurate  when  these  boundaries  are  close  to  the  target,  but 
not  close  enough  to  be  included  in  the  computational  domain. 

5.  Due  its  modular  nature,  the  IMP  technique  lends  itself  to  the  use  of  parallel  processing, 
including  GPU  processing  and  fast  multipole  method  (FMM)  [14]-[18]  to  substantially  im¬ 
prove  its  speed  even  further. 


Results  and  Discussions 

In  this  section,  we  validate  the  axially-symmetric  and  3D  versions  of  our  method,  Axi-IMP  and 
3D-IMP  using  analytical  solutions  when  possible,  and  semi-analytic  solutions  based  on  the  T- 
Matrix  method.  The  use  of  analytic  solutions  as  benchmark  solutions  is  motivated  by  the  desire 
to  demonstrate  the  accuracy  of  our  method  without  having  to  worry  about  the  accuracy  of  the 
benchmark  solution.  Obviously,  analytical  solutions  are  numerically  exact  and  even  though  the 
T-Matrix  solution  in  only  semi-analytic,  it  is  the  next  best  solution,  particularly  for  a  sphere.  After 
the  models  are  validated,  we  apply  them  to  compute  the  acoustic  color  for  various  UXOs  and 
compare  results  with  measurements. 
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Applications 

Axi-IMP 

We  applied  the  axially  symmetric  formulation  (Axi-IMP)  to  compute  the  acoustic  color  templates 
for  a  variety  of  targets  and  compared  the  execution  times  between  this  technique,  the  3D-IMP 
and  the  axially  symmetric  finite  element  solution,  AxiScat.  Figure  (10)  shows  the  results  for  an 
aluminum  replica  of  a  UXO,  whose  length  is  approximately  48  cm  and  its  maximum  diameter  is 
about  10  cm.  The  aspect  angle  is  measured  from  the  tip  the  UXO,  where  angles  0,  90  and  180 
degrees  correspond  to  tip-on,  broadside  and  end-on  incidence,  respectively.  The  top  left  panel  was 
obtained  by  the  axially  symmetric  finite  element  model,  AxiScat,  and  it  took  40  hours.  The  top 
right  figure  was  obtained  by  the  Axi-IMP  and  it  only  took  3  hours  for  10  azimuthal  orders.  The 
same  computation  was  carried  out  by  3D-IMP  (not  shown)  and  that  took  24  hours.  There  is  a 
rule  of  thumb  for  determining  the  minimum  number  of  azimuthal  orders  needed  for  a  solution  to 
converge.  For  targets  similar  in  size  and  the  range  of  frequencies  shown  here,  this  number  is  on 
the  order  of  20.  However,  we  used  azimuthal  orders  that  were  less  than  20  and  the  solutions  seem 
to  have  converged. 


Figure  10:  The  acoustic  color  template  for  the  aluminum  UXO  replica.  The  top  left  panel  was 
obtained  using  AxiScat  and  it  took  40  hours.  The  top  right  panel  was  obtained  using  Axi-IMP  and 
it  took  only  3  hours  for  10  azimuthal  orders.  The  bottom  left  panel  shows  the  2D  quadrilateral 
mesh  on  the  cross  section  of  the  UXO,  which  is  revolved  around  the  axis  of  symmetry  to  obtain  the 
3D  object.  The  bottom  right  panel  is  a  picture  of  the  target. 


The  acoustic  color  templates  for  another  UXO  referred  to  as  Bullet  105  are  shown  in  Fig.  (11). 
The  aspect  angle  in  this  simulation  is  measured  from  broadside,  where  -90  degrees  corresponds 
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to  tip-on  incidence  and  90  degrees  corresponds  to  end-on  incidence.  The  length  of  this  target  is 
slightly  under  50  cm  and  its  maximum  exterior  diameter  is  approximately  1 1  cm.  It  is  made  of 
three  materials,  aluminum,  steel  and  copper.  In  the  bottom  right  panel  of  Fig.  (11)  the  shape  of 
the  Bullet  105  is  depicted,  where  red,  gray  and  light  red  represents  copper,  steel  and  aluminum, 
respectively.  It  also  has  a  cavity  inside  that  in  this  simulation  was  assumed  to  be  filled  with  air.  The 
axially  symmetric  finite  element  computation  for  this  target  using  AxiScat,  shown  in  the  top  left 
panel,  was  carried  out  by  APL-UW,  so  we  are  not  sure  how  long  it  took,  but  our  computation  using 
the  Axi-IMP,  shown  in  the  top  right  panel,  only  took  3.8  hours  for  16  azimuthal  orders  and  based 
on  the  results  shown  in  Fig.  (10),  we  can  safely  assume  that  our  solution  is  at  least  10  times  faster  2. 


105  APL-UW  Solution 
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Figure  11:  The  acoustic  color  template  for  Bullet-105.  The  top  left  panel  was  obtained  using  Ax¬ 
iScat  .  The  top  right  panel  was  obtained  using  AxiSym- IMP  and  it  took  3.8  hours  for  16  azimuthal 
orders.  The  bottom  left  panel  shows  the  2D  quadrilateral  mesh  on  the  cross  section  of  the  target, 
which  is  revolved  around  the  axis  of  symmetry  to  obtain  the  3D  object.  The  bottom  right  panel 
shows  a  depiction  of  the  target. 


We  next  computed  the  acoustic  color  template  for  the  Howitzer  UXO.  This  is  a  complex  target 
made  of  five  materials,  steel,  mild  steel,  copper,  plastic  and  aluminum.  This  target  is  about  85  cm 
long  and  its  maximum  radius  is  about  10  cm.  In  our  simulations  we  assumed  that  the  target  was 
filled  with  air.  The  acoustic  color  template  took  9.5  hours  to  compute  using  16  azimuthal  orders. 
The  results  are  shown  in  Fig.  (12). 

2We  compared  the  execution  times  between  AxiScat  and  Axi-IMP  by  computing  the  acoustic  color  for  a  1  x  2  foot 
aluminum  cylinder  from  0  to  10  kHz  and  from  0  to  90  degrees  on  a  Mac  Pro.  AxiScat  took  8  hours,  while  Axi-IMP 
took  only  19  minutes. 
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Figure  12:  The  acoustic  color  template  for  the  Howitzer  UXO.  The  tip-on  incidence  is  at  -90 
degrees  and  the  end-on  incidence  is  at  90  degrees.  The  figure  on  the  far  right  shows  the  finite 
elemen  t  mesh  used  in  this  computation.  The  computation  took  9.5  hours  for  16  azimuthal  orders. 

The  above  run  times  for  the  models  used  are  summarized  in  the  table  below.  Note  that  ’N/A’  in  the 
table  indicates  that  we  have  no  data  on  the  corresponding  models/targets,  but  based  on  the  other 
results,  rough  estimates  of  run  times  can  be  assessed.  Note  also  that  even  though  solution  in  free 
space  is  listed  as  a  limitation  for  both  AxiScat  and  Axi-IMP,  solution  in  a  waveguide  can  be  ob¬ 
tained  by  adding  multi-paths  in  post  processing.  This  is  done  in  the  AxiScat  examples  shown  in  this 
report  (see  Figs.  14-16)  and  can  also  be  done  using  the  Axi-IMP  solution.  But  these  solutions  are 
approximate  and  not  self-consistent  since  the  finite  element  solution  for  each  multi-path  assumes 
that  the  target  is  in  free  space.  The  3D-IMP  solution  does  not  have  any  of  these  limitations. 


AxiScat 

Axi-IMP 

3D-IMP 

Aluminum  Cylinder 

8  hours 

19  minutes 

1.5  hours 

Aluminum  UXO 

40  hours 

3  hours 

24  hours 

Bullet-105 

N/A 

3.8  hours 

32  hours 

Howitzer  UXO 

N/A 

9.5  hours 

85  hours 

Limitations 

Axi-symmetric/Free  Space 

Axi-symmetric/Free  Space 

None 

Table  1:  Comparison  of  run  times  between  different  models.  Note  that  the  run  times  for  the  alu¬ 
minum  cylinder  are  for  100  frequencies  between  1-10  kHz.  For  the  other  targets  they  are  for  300 
frequencies  between  1-30  kHz.  Also  note  that  the  run  times  for  the  3D-1MP  model  is  competitive 
with  the  other  two  models  despite  the  fact  that  it  does  not  take  advantage  of  the  axially -symmetric 
nature  of  the  target,  as  do  the  other  two  models.  ’N/A’  indicates  that  no  data  are  available  for  the 
corresponding  models/targets. 
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3D-IMP 


We  pointed  out  earlier  that  despite  its  amazing  speed,  the  Axi-IMP  model  can  only  be  applied  to 
axially- symmetric  targets  in  free  space.  For  general- shaped  targets  or  even  for  axially-symmetric 
targets  in  an  environment  other  than  free  space,  we  must  use  the  3D  version  of  the  model,  3D-IMP. 
In  what  follows  we  apply  this  model  to  simulate  various  experimental  scenarios.  One  such  scenario 
is  depicted  by  Fig.  (13),  which  describes  the  experimental  setup  involving  an  aluminum  replica  of 
a  proud  UXO  used  in  one  of  the  PondEx  experiments  off  the  coast  of  Panama  City,  Florida. 


Figure  13:  The  PondEx  experiment  involving  scattering  from  a  proud  UXO. 

In  this  experiment,  the  target  is  proud  on  the  bottom  and  the  source  and  receiver  trace  a  semi-circle 
of  4.5  m  radius  around  the  target  at  a  plane  3.6  m  from  the  bottom.  In  this  setup  zero  and  180 
degrees  corresponded  to  tip-on  and  end-on  incidence.  In  our  simulation,  we  modeled  the  ocean 
environment  as  two  fluid  half  spaces.  The  model  values  for  the  bottom  density  and  sound  speed 
were  2000  kg/m3  and  1694  m/s,  respectively.  The  acoustic  color  is  shown  in  Fig.  (14),  where 
the  left  panel  shows  the  experimental  results,  the  middle  panel  our  results  and  the  right  panel 
results  obtained  by  APL-UW.  We  see  that  the  two  modeled  results  agree  very  well  and  they  agree 
reasonably  well  with  the  experimental  results. 


0  10  20  30  10  20  30  10  20  30 

Frequency  (kHz)  Frequency  (kHz)  Frequency  (kHz) 


Figure  14:  The  acoustic  color  results  for  the  proud  aluminum  UXO.  The  left  panel  shows  exper¬ 
imental  results,  the  middle  computed  using  our  3D  model,  3D-IMP,  and  the  right  panel  shows 
computed  results  from  APL-UW  using  their  Hybrid  method. 
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Similar  results  for  a  buried  UXO  with  otherwise  identical  experimental  geometry  are  shown  in 
Fig.  (15).  The  target  is  buried  in  such  a  way  that  its  axis  of  symmetry  is  parallel  to  the  water/bottom 
interface  and  its  top  is  approximately  5  cm  below  the  interface.  Again  there  is  good  agreement 
between  the  two  models  and  reasonable  agreement  with  the  measured  results. 
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Figure  15:  The  acoustic  color  results  for  the  buried  aluminum  UXO.  The  left  panel  shows  exper¬ 
imental  results,  the  middle  computed  using  our  3D  model,  3D-IMP,  and  the  right  panel  shows 
computed  results  from  APL-UW  using  their  Hybrid  method. 


The  acoustic  color  for  a  half  buried  UXO  is  shown  in  Fig.  (16),  where  in  this  case  the  source/receiver 
geometry  is  the  same,  except  the  UXO  is  buried  in  such  a  way  that  its  axis  of  symmetry  is  parallel 
to  and  at  the  same  level  as  the  interface.  Unfortunately,  we  were  not  able  to  obtain  measured  results 
for  this  case,  so  we  compare  our  modeled  results  with  those  of  APL-UW’s  Hybrid  technique.  It 
should  be  pointed  out  that  the  small  differences  between  our  results  and  those  of  APL-UW  may  be 
because  our  solution  is  numerically  exact  and  thus  includes  multiple  scattering  to  all  orders,  where 
those  of  APL-UW  accounts  for  only  four  principal  rays. 

In  examining  the  acoustic  color  for  the  aluminum  replica  of  a  UXO  shown  in  Fig.  (10)  we  notice 
the  presence  of  three  resonances  around  3,  6  and  9  kHz  and  aspect  angles  of  40  and  140  degrees. 
These  resonances  are  due  to  the  first  three  bending  modes  of  the  UXO.  The  first  bending  mode  is 
a  high  Q  mode  that  occurs  around  3  kHz.  The  reason  these  resonances  occur  at  the  above  angles 
is  that  the  UXO  radiates  most  efficiently  at  these  angles.  Simulations  show  that  this  seems  to  be 
the  characteristic  of  all  objects  of  this  size  with  a  circular  cross  section.  Examining  the  acoustic 
color  for  the  same  UXO  in  a  two  half-space  environment  in  Figs.  (14),  (15)  and  (16)  shows  the 
existence  of  the  same  resonances,  but  with  angles  shifted  towards  tip-on  and  endfire  directions  by 
the  presence  of  the  interface.  However,  these  resonances  seem  to  be  robust  features  that  can  be 
used  for  classification. 
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Figure  16:  The  acoustic  color  results  for  the  half  buried  aluminum  UXO.  The  left  panel  shows 
computed  results  using  our  3D  model,  3D-IMP,  and  the  right  panel  shows  those  from  APL-UW 
using  their  Hybrid  method. 

Since  the  computation  of  matrices  A  and  B  for  the  exact  Green’s  functions  given  by  (31)  and  (32) 
is  numerically  very  intensive,  it  has  been  common  practice  to  approximate  these  Green’s  functions 
by  free  space  Green’s  functions.  This  is  called  the  single-scatter  approximation.  We  pointed  out 
before  that  if  the  Green’s  functions  do  not  satisfy  the  waveguide  boundary  conditions,  an  exact 
solution  cannot  be  obtained,  and  among  other  things,  the  solution  will  not  be  able  to  account  for 
multiple  scattering.  To  experimentally  test  the  accuracy  of  the  single-scatter  approximation,  an 
experiment  was  carried  out  at  Washington  State  University. 

In  the  experiment  the  details  of  which  are  given  in  [19],  an  aluminum  cylinder  ( length  —  50.8  mm, 
radius  =  12.7  mm)  was  suspended  from  two  strings  (not  shown)  attached  to  the  top  and  bottom 
of  the  cylinder.  The  lengths  the  strings  could  vary,  thereby  allowing  the  tilt  angle  of  the  cylinder 
to  vary.  The  strings  were  attached  to  a  circular  platform  above  the  water,  which  could  be  turned, 
thus  allowing  the  cylinder  to  be  rotated  in  a  horizontal  plane.  The  transducer  was  stationary  in 
the  water,  producing  an  incident  field  that  made  a  19.5°  angle  with  the  water  surface.  By  rotating 
the  cylinder,  the  incident  angle  could  be  varied,  exposing  the  cylinder  to  the  incident  field  in  a  full 
angular  range  of  180°.  In  the  scenario  shown  in  the  figure,  the  lengths  of  the  strings  are  adjusted 
in  such  a  way  that  the  side  attached  to  one  string  is  exactly  at  air-water  interface  and  the  axis  of  the 
cylinder  makes  a  34°  angle  with  the  surface  of  the  water.  In  our  simulations,  instead  of  rotating  the 
cylinder,  we  rotated  the  source  as  is  shown  in  Fig.  (17). 
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Figure  17:  The  tilted  cylinder  experiment  geometry. 

In  this  geometry,  when  the  source  is  at  6  —  90°  the  incident  field  is  multiply  scattered  by  the 
cylinder  and  the  water  surface.  This  is  shown  in  Fig.  (18) 


Figure  18:  Multiple  scattering  by  the  cylinder  and  the  water  surface. 

The  measured  acoustic  color  plot  for  the  tilted  cylinder  is  shown  in  Fig.  (19),  where  the  strong 
response  due  to  the  multiply-scattered  ray  can  be  seen  at  0  —  90°  for  ka  f  6.5. 
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Figure  19:  The  measured  acoustic  color  as  a  function  ofka  and  angle  0  (see  Fig.  (17)).  The  strong 
response  due  to  the  multiply -scattered  ray  can  be  seen  at  9  —  90°  and  ka  6.5. 

Since  the  pressure  is  zero  at  the  water  surface  (z  =  0),  the  Green’s  function  satisfying  the  boundary 
condition  at  the  surface  can  be  obtained  from  the  method  of  images 

G{x,y,z;,x' ,y  ,z')  =  ^{G  -G+},  (34) 
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where 


and 


G  (x,y,z;x' ,y' ,z') 


eikyJ  ( x-x')2+(y~y)2+{z-z ')2 

sj (x  -  x')2  +  (y  -  y')2  +  (z  —  z')2’ 


G+  (x,y,z;,x',y',zr) 


eikyJ  (.v-y)2+(.v-y')2+(z+z')2 
y/(x-.r/)2  +  (y-y)2  +  (z  +  z')2 


The  acoustic  color  computed  using  our  model,  which  uses  the  correct  Green’s  function  given  by 
(34)  and  the  single-scatter  model,  which  uses  a  free  space  Green’s  function  (the  first  term  in  (34)) 
are  shown  in  Fig.  (20). 


Tilt=34°  (3D-IMP)  Tilt=34°  (Single-Scatter  Approx.) 


Figure  20:  The  acoustic  color  for  the  tilted  cylinder  computed  using  our  3D-IMP  model  on  the  left 
and  a  model  that  uses  the  single- scatter  approximation  on  the  right.  The  strong  response  due  to 
the  multiply -scattered  ray  can  be  seen  at  0  =  90°  and  ka  f  6.5  in  the  3D-IMP  solution,  but  it  is 
absent  in  the  single-scatter  solution. 


This  figure  shows  that  the  strong  response  due  to  the  multiply-scattered  ray  at  0  —  90°  is  present 
in  the  result  on  the  left,  which  uses  the  correct  Green’s  function,  and  is  absent  in  the  single-scatter 
solution.  This  proves  that  to  be  able  account  for  multiple  scattering,  the  single-scatter  approxima¬ 
tion  is  not  adequate  and  one  should  use  the  proper  Green’s  functions,  i.e.  those  that  satisfy  the 
environment  boundary  conditions. 

As  a  final  example  that  involves  the  use  of  both  Axi-IMP  and  3D-IMP  models,  we  computed 
the  acoustic  color  for  the  Bullet-105  UXO  with  a  hole  drilled  onto  its  side.  This  allows  water  to 
enter  the  UXO,  which  will  affect  its  acoustic  response.  Further,  from  a  computational  point  of 
view,  the  presence  of  the  hole  breaks  the  axial  symmetry  of  the  target  and  the  problem  can  no 
longer  be  solved  by  an  axially- symmetric  model.  So  we  used  our  3D  model,  3D-IMP.  With  the 
computational  advantages  that  3D-IMP  has,  it  was  able  to  solve  this  problem  in  70  hours  on  an 
Apple  Mac  Pro.  With  a  conventional  finite  element  method  it  is  impossible  to  solve  this  problem 
on  a  regular  desktop  computer.  The  results  are  shown  in  Fig.  (21).  The  figure  on  the  left  shows  the 
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3D-IMP  solution  for  the  UXO  with  a  hole,  and  for  comparison,  the  figure  on  the  right  shows  the 
acoustic  color  for  the  same  UXO  without  the  hole,  which  we  obtained  using  our  axially-symmetric 
model,  Axi-IMP.  This  solution  only  took  3.8  hours  to  run,  which  shows  once  more  the  compu¬ 
tational  advantage  of  the  axially-symmetric  model,  when  it  can  be  employed.  The  bottom  two 
figures  show  CAD  illustrations  of  the  UXO  with  and  without  the  hole.  It  is  important  to  note  that 
the  presence  of  the  hole  changes  the  acoustic  color  significantly. 


Figure  21:  The  acoustic  color  for  a  the  Bullet-105  UXO  with  a  hole  drilled  on  its  side  on  the 
left  and  the  same  solution  for  the  UXO  without  the  hole.  The  angle  0  corresponds  to  tip-on,  90 
to  broadside  and  180  to  end-on  incidence.  The  bottom  two  figures  show  CAD  illustrations  of  the 
UXO  with  and  without  the  hole. 

Summary  and  Conclusions 

During  the  last  three  years,  we  developed  numerical  tools  to  compute  scattering  from  an  arbitrary 
elastic  target  in  an  arbitrary  ocean  environment.  These  models  are  based  on  the  fluid-structure 
interaction  method,  which  uses  the  finite  element  method  to  compute  the  impedance  matrix  for  the 
target  and  the  boundary  element  method  to  compute  the  propagation  in  the  environment  in  which 
the  target  is  embedded.  The  two  solutions  are  coupled  by  imposing  the  continuity  of  the  pressure 
and  normal  particle  velocity  on  the  surface  of  the  target.  This  method  produces  a  system  of  equa¬ 
tions  that  is  self-consistent  and  rigorously  integrates  propagation  and  scattering.  This  means  that 
it  produces  a  numerically  exact  solution  for  scattering  from  an  elastic  target  in  any  ocean  environ¬ 
ment  that  can  be  characterized  by  a  Green’s  function. 

We  also  developed  an  axially-symmetric  version  of  the  model,  which  we  refer  to  as  the  Axi-IMP. 
This  model  can  be  applied  to  an  axially-symmetric  target  for  an  arbitrary  (not  necessarily  axially- 
symmetric)  incident  field  in  free  space.  Since  most  of  the  targets  of  interest  are  axially-symmetric, 
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this  model  proves  to  be  very  useful.  It  was  shown  in  this  report  that  scattering  from  an  axially- 
symmetric  target  in  a  two  half-space  environment  can  be  solved  by  using  a  combination  of  ray 
theory  to  determine  which  rays,  in  addition  to  the  incident  ray,  are  incident  on  the  target,  and  the 
AxiScat  model  to  compute  the  scattered  field  for  each  ray.  The  problem  can  then  be  solved  as 
a  superposition  of  the  solution  for  each  ray.  Even  though  this  is  not  an  exact  solution,  it  is  very 
accurate,  as  was  shown  in  Figs.  (14),  (15)  and  (16).  Our  axially-symmetric  model,  Axi-IMP,  can 
do  these  types  of  computations  at  much  faster  speeds  (see  the  footnote  on  page  21). 

We  validated  both  models  by  computing  the  backscattered  field  from  a  sphere  in  free  space  and 
compared  our  results  with  those  obtained  from  the  classical  partial  wave  solution  of  the  sphere.  We 
also  validated  the  3D-IMP  by  computing  backscattering  from  a  proud,  partially  buried  and  fully 
buried  sphere  and  compared  our  results  with  those  of  the  T-Matrix  method.  Once  the  models  were 
validated,  we  computed  the  acoustic  color  for  various  UXOs  in  free  space  and  from  an  aluminum 
replica  of  a  UXO  that  was  used  in  the  PondEx  experiments.  In  the  latter  case,  we  compared  our 
results  to  measurements  and  showed  that  the  our  models  matched  experimental  results  pretty  well 
and  contained  most  of  the  features  seen  in  the  measured  results.  We  also  demonstrated  that  the  30- 
IMP  correctly  accounted  for  multiple  scattering,  as  advertised,  by  computing  the  acoustic  color  for 
a  1  x  2  inch  aluminum  cylinder  near  the  air-water  interface  and  comparing  our  results  with  mea¬ 
surement.  These  measurements  were  made  in  a  carefully  designed  experiment  at  Washington  State 
University  to  look  for  the  effects  of  multiple  scattering.  Finally,  we  applied  the  3D-IMP  to  compute 
scattering  from  a  UXO  with  a  hole  drilled  on  its  side.  Even  though  the  UXO  is  axially-symmetric, 
the  presence  of  the  hole  breaks  this  symmetry  and  the  problem  can  only  be  solved  by  a  3D  model. 
We  compared  the  acoustic  color  for  this  case  with  that  without  the  hole  and  showed  that  the  pres¬ 
ence  of  the  hole  significantly  changes  the  acoustic  color,  not  particularly  because  of  the  presence 
of  the  hole,  but  because  the  UXO  is  flooded.  This  is  the  type  of  simulation  that  would  be  needed 
to  study  the  effects  of  long-term  exposure  to  the  harsh  ocean  environment  on  the  acoustic  response 
of  the  UXOs. 

Even  though  the  3D-IMP  can  compute  the  acoustic  color  for  all  UXOs  of  interest  in  reasonable 
times,  the  execution  times  can  be  reduced  drastically  by  the  use  of  modern  techniques  like  the  fast 
multipole  method  (FMM)  [14]-[18]  and  hierarchical  matrices  (H-Matrices)  [20]  and  [21].  One 
reason  that  the  current  execution  times  are  manageable  is  that  problems  are  meshed  at  the  lowest 
levels  possible.  For  example,  the  rule  of  thumb  is  that  one  should  have  on  the  order  15  elements 
per  wavelength  for  an  accurate  solution,  but  lack  of  memory  and  CPU  force  this  number  to  be  on 
the  order  6-10.  The  use  of  the  above-mentioned  techniques  reduces  most  operations  that  are  on  the 
order  0(n2)  to  0(n  x  log(n)),  where  n  is  the  number  of  elements.  This  is  achieved  by  rigorously 
maintaining  a  pre-specified  error  tolerance.  With  type  of  reduction  in  execution,  one  can  be  more 
generous  with  the  number  of  elements  and  thus  also  increase  the  accuracy  of  the  solution.  We  plan 
to  implement  these  techniques  in  our  models  in  the  future. 
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